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Abstract 

We present a non-perturbative technique to study pulse dynamics in excitable me- 
dia. The method is used to study propagation failure in one-dimensional and two- 
dimensional excitable media. In one-dimensional media we describe the behaviour 
of pulses and wave trains near the saddle node bifurcation, where propagation fails. 
The generalization of our method to two dimensions captures the point where a 
broken front (or finger) starts to retract. We obtain approximate expressions for 
the pulse shape, pulse velocity and scaling behavior. The results are compared with 
numerical simulations and show good agreement. 
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Excitable media are often found in biological and chemical systems. Ex- 
amples of excitable media include electrical waves in cardiac and nerval tissue 
PQ, [2], cAMP waves in slime mold aggregation [3j and intracellular calcium 
waves pi]. Excitable media support localized pulses and periodic wave trains. 
In 2 dimensions rotating vortices (or spirals) are possible [5j. The critical 
behavior of pulses, wave trains and spirals, i.e. propagation failure, is often 
associated with clinical situations. The study of spiral waves is particularly 
important as they are believed to be responsible for pathological cardiac ar- 
rhythmias Spiral waves may be created in the heart through inhomo- 
geneities of the properties of the cardiac tissue. 

We investigate critical behavior relating to these 3 wave types. We develop a 
non-perturbative test function method which allows to study the bifurcation 
behavior of critical waves. In particular, we study under what conditions a 
broken front will sprout and develop into a spiral wave or retract. Analyti- 
cal formulas for the growing velocity of a broken front are given. For wave 
trains we provide a time dependent extension which supports a Hopf bifurca- 
tion which is also observed in numerical simulations of excitable media. This 
seems to be related to alternans OIHl, which also are discussed in the context 
of cardiac electric pulse propagation. The methods and results are general, 
and can be applied to other excitable media. 

1 Introduction 

Many chemical and biological systems exhibit excitability. In small (zero-dimensional) ge- 
ometry they show threshold behavior, i.e. small perturbations immediately decay, whereas 
sufficiently large perturbations decay only after a large excursion. One dimensional (ID) 
excitable media support travelling pulses, or rather, periodic wave trains ranging in wave- 
length L from the localized limit L ^ oo to a minimal value L^.. Pulses and wavetrains 
are best-known from nerve propagation along axons. In 2D one typically observes spiral 
waves. Spirals have been observed for example in the auto-catalytic Belousov-Zhabotinsky 
reaction [S] , in the aggregation of the slime mold dictyostelium discoideum [3] and in car- 
diac tissue 12]. 

For certain system parameters the propagation of pulses, wave trains or the develop- 
ment of spiral waves may fail (see for example [HE]). The analytical tools employed to 
describe these phenomena range from kinematic theory fjl^, asymptotic perturbation 
theory to dynamical systems approaches [TH] I17j. However, no theory exists 

which describes propagation failure using only equation parameters, and which repro- 
duces the behavior close to the bifurcation point. For example, asymptotic perturbation 
theory fails to describe the square-root scaling behavior of the amplitude and the pulse 
velocity with respect to the bifurcation parameter at the bifurcation point. In kinematic 
theory results are not given entirely in terms of the system parameters. In this paper we 
develop a non-perturbative method to study propagation failure and compare the results 
with numerical simulations. 
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Most theoretical investigations are based on coupled reaction-diffusion models. We 
follow this tradition and investigate a two- component, two-dimensional excitable medium 
with an activator u and a non-diffusive inhibitor v described by 

dfU = DAu + J^{u,v), J^{u,v) = u{l — u){u — Us — v) 

dtv = e {u — a v) . (1) 

This is a reparametrized version of a model introduced by Barkley ^Hl- We expect our 
method to be independent of the particular model used. Note that the diffusion constant 
D is not a relevant parameter as it can be scaled out by rescaling the length. 

This model incorporates the ingredients of an excitable system in a compact and lucid 
way. Thus, for > the rest state Mq = = is linearly stable with decay rates cxi = Ug 
along the activator direction and a2 = ea along the inhibitor direction. Perturbing u above 
the threshold Ug (in OD) will lead to growth of u. In the absence of v the activator would 
saturate at m = 1 leading to a bistable system. A positive inhibitor growth factor e and 
a > forces the activator to decay back to m = 0. Finally also the inhibitor with the 
refractory time constant (e a)~^ will decay back to u = 0. For a > 1/(1 — Ug) the system 
is in zero-dimensional systems no longer excitable but instead bistable. This relaxation 
behaviour in the OD system for super-threshold perturbations gives rise to pulse solutions 
in the ID (and 2D) case. The relaxation mechanism mediated by the inhibitor forces the 
pulse solution to decay in its back. Hence we observe pulses in excitable media and not 
fronts. 

The Barkley model is a variant of the class of 2-component Fitzhugh-Nagumo models. 
In the traditional Fitzhugh-Nagumo models J^{u,v) in Eq.((TI) is replaced by J^p]\r{u,v) = 
u{l — u){u — Us) — V. Thus the nuUclines T{u^ v) = 0, which in the traditional Fitzhugh- 
Nagumo models are cubic polynomials, are replaced by straight lines in the Barkley model. 
Most of the qualitative behavior in the relevant parameter ranges is unchanged by this. 
Whereas the Barkley model is computationally more efficient and also analytically better 
tractable, the traditional Fitzhugh-Nagumo models display a feature which makes them 
more realistic for the discription of excitable media in biology: the activator experiences 
an undershoot below its equilibrium value and slow decay in the tail region of a pulse. 
We expect that for the phenomena discussed in this paper the difference only leads to 
quantitative changes (indeed, we have done some tests to verify this assertion). 

In order to study pulse propagation in ID it is useful to first consider the case of 
constant v. The resulting bistable model is exactly solvable fH] and the pulse velocity is 

Cf{v) = \J^[^ — '^{us + v)]. Hence, excitability requires that Ug is below the stall value 

|. The quantity IS. = \ — Ug characterizes the strength of excitability and c/(0) coincides 
with the solitary pulse velocity for e — * 0. 

Clearly, for Ug < Uc = \ and not too large a, pulse propagation fails for e larger than 
some ec- The critical growth factor ec describes the onset of a saddle-node bifurcation 
^UE]. The saddle node can be intuitively understood when we consider the activator 
pulse as a heat source, not unlike a fire-front in a bushfire. Due to the inhibitor the 
width of the pulse decreases with increasing e. Hence, the heat contained within the pulse 
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decreases. At a critical width, or a critical e, the heat contained within the pulse is too 
small to ignite/excite the medium in front of the pulse. 

For periodic wavetrains the saddle node depends on the wavelength. The pulses run 
into the inhibitor field of their respective preceding pulse. Hence, propagation failure for 
periodic wave trains is controlled by the decay of the inhibitor and propagation is only 
possible when the inter-pulse distance L becomes larger than a critical wavelength L^. 
Note that diverges for a ^ when the decay rate of the inhibitor cr2 vanishes. 

In previous analytic works on one-dimensional pulses the limit of small e was considered 
PT| IT^ . Then the solution of the activator u is well separated in two fiat plateau regions 
with u ^ 1 and m = which are separated by a steep narrow front. This approach does 
not capture the square- root behavior of co(e) close to the saddle-node bifurcation point. 
Close to the bifurcation point the solution resembles rather a bell-shaped pulse than a 
plateau. In this paper we will be concerned with the behavior near the bifurcation at ec 
and make explicit use of the observed bell-shape form of the solution. This allows us to 
describe the scaling behavior close to the bifurcation in terms of the equation parameters. 

In 2 dimensions spiral waves are observed. Spirals can be created in excitable media 
from a finger, i.e. a ID pulse which is extended in the second dimension and has one free 
end. Fingers may be created due to inhomogeneities in the excitability of the system |2Uj . 
At long times, the free end will either sprout or retract depending on the growing velocity 
Cg being positive or negative. If Cg > 0, the tip of the finger will sprout into the fresh 
medium, and in particular it sprouts and curves backwards. This causes a non-vanishing 
curvature at the tip of the finger. Due to the increased curvature the finger tip is slower 
than the fiat part of the front further away from the tip. Thus, the extending part far 
from the tip will curl up. This leads to the formation of a spiral with the free end at its 
core. The criterion > is therefore a necessary criterion for spiral formation. 

The transition to retraction always occurs before the ID propagation failure. It is 
harder to tackle analytically. General and universal dynamical system approaches ex- 
ploiting only the Euclidean symmetry of excitable media describe the transition as a drift 
bifurcation. These model independent theories give an explanation for the divergence of 
the core radius at the bifurcation point and also explain why a finger at the bifurcation 
point is translating with finite speed, i.e. that the transition to retracting fingers always 
occurs before the ID propagation failure ^7j. Unfortunately, it has the unphysical re- 
sult that at the bifurcation point eg a spiral changes its sense of rotation ^Tj. Whereas 
these theories treat the spiral as a global solution of the underlying equations and see 
the transition to spiral waves as a pitchfork or drift bifurcation [THl I17j . we take a local 
approach and describe the transition not as a bifurcation but as a quantitative change in 
the velocity of the finger, analogous to a Maxwell point in a first-oder phase transition. 
Asymptotic techniques in the limit e <C 1 have also been performed for this problem 
[T^ [T3] and produced an analytical expression for the onset of retraction for small e. 
Moreover, the authors were able to go one step further, in a detailed and sophisticated 
asymptotic analysis, and described the onset of meandering. In this paper we go beyond 
the restriction of small e and propose a more intuitive approach to the problem of the 
growing velocity which nevertheless gives good agreement with the numerics. 
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In the following we assume given excitability parameters a, Us and take e as the 
bifurcation parameter. We will be looking for solutions that move with velocity Cq in the 
X direction and may grow (or retract) with velocity Cg in the y direction. Thus we rewrite 
Eqs. (0) in a frame moving with velocity (cq, Cg) as 

D{dl + dl)u + cod,u + CgdyU + J^{u,v) = (2) 
CodxV + CgdyV + e{u — a v) = . (3) 



2 Pulse propagation in one-dimensional excitable me- 
dia 

We first look for pulse and wavetrain solutions that do not depend on y. The reason 
for the failure to describe the pulse properties at ec within the framework of asymptotics 
employing the smallness of e is due to the fact that at ec the pulse shape for u cannot be 
separated into a steep narrow front and a flat plateau. Hence, asymptotic techniques such 
as inner and outer expansions are bound to fail. Instead the pulse has the shape of a rather 
symmetric bell-shaped function (see Fig. 1). In the following we make explicit use of the 
shape of the pulse close to the critical point and parametrize the pulse appropriately; a 
method reminiscent of the method of collective coordinates in the studies of solitary waves 

m 

We choose u of the general form 

u{x) = foU{r]) with r] = wx , (4) 

where U{ri) is chosen as a symmetric, bell-shaped function, for example a Gaussian, of 
unit width and height. Hence, we restrict the solutions to a subspace of bell-shaped 
functions U (77) which is parametrized by the amplitude /o and the inverse pulse width w. 
The aim of our method is to determine the so far undetermined parameters. This is done 
by minimizing the error made by the restriction to the subspace defined by (jH). 

We avoid further uncontrolled approximations and solve for the inhibitor field v ex- 
plicitly 



;(ry) = foGVir]) with V{r]) = e"^" (v* - j^ri'e-"^"' U{ri')^ , 



(5) 



where 

e = e/{cow) . (6) 
y* is determined via the periodic boundary condition v{—wL/2) = v{wL/2) and is 

1 /"■§■"' 

V' = , . , / dv'e'^'^'-^'-^-^Uiv') . (7) 

2smh(aB^u;) J-L^ 
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We assume that the width is small compared to the distance between two consecutive 
pulses L. in the temporal domain this means that the time scale for the decay of the 
inhibitor is much longer than the activator pulse width. This assures that the activator 
field of consecutive pulses is well separated and only the inhibitor field overlaps, and 
that the interaction between pulses is mediated only through the inhibitor. Otherwise we 
would have to choose periodic functions U{ri). However, in this case we can now replace 
the limits of integration by ±00. For the isolated pulse, i.e. when L 00, V* 

vanishes. 

We now determine the parameters /o and w by projecting Eq. Q onto the tangent 
space of the restricted subspace defined by (jH). The tangent space is spanned by du/dfo = 
U and du/dw = rjUr). This assures that the error made by restricting the solution space 
to the test functions is minimized. To achieve this, we multiply Eq. ((2)) with the basis 
functions of the tangent-space U and rjUj^, integrate over the //-domain and require the 
projection to vanish, i.e. 

{Dw'^Unn + u{l -u){u-Us- v)\U)u=foU(ri) = (8) 
{Dw'^Unrj + u{l -u){u-Us- v)\vU')u=foU{-n) = , (9) 

where the brackets indicate integration over the whole ?7-domain. The terms proportional 
to the velocity cq vanish. 

The resulting equations can be combined to give, at fixed and a, a quadratic equation 
for /o with two solutions /o± which describe the stable and unstable branch, respectively 
(see below for a subtle issue at the saddle node). We obtain 

Af^ + Bfo + C = 0, (10) 

where 

A = |(f/4) - f (f/3y) - ^{7]U^V) , 
B = -|(1 + n,)(f/3) + e(f/V) + ^{rjU^V) , C = Us{U^) . (11) 

The corresponding inverse width parameters w± for the stable and unstable branch are 
given by 



w 



[f^{-{u') + e(f/V)) + /o((i + us){u') - e(f/V)) - us{u')] . (12) 



The velocity cq can now be determined in the standard way by multiplying Eq. ((2|) by 
and integrating over x. We obtain 

poo 

/ ^ut.dx 

J —00 



/o© r/o//rr4\ ^o/r^3TA^ ^^/rrSx ^o/rr2i 



w{U^) 3 2 

Finally, we can determine e from e = cqwQ. Multiplying (fT!^ by w one sees that e can be 
computed without calculating w and cq. 
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We will use a Gaussian U = e"'' as an ansatz function Note that other symmetric 
bell-shaped functions such as a sec/i-function, are possible, too. Then one has (U^) = 
-y/vr/n and (U^) = (f/^) = \/tt/2. The parameters of the test function /o and w, and the 
front velocity cq can now be determined numerically using ()10I12I13|) . 

Simplifications are possible in the useful limit Qa ^1. In this limit the (temporal) 
inverse pulse width (ty„c„)~^ is small compared to the inhibitor decay time (ea)~^ (see 
definition (jH)). Then in Eqs. and (jl3|) the terms proportional to a can be omitted, 
and for the calculation of {If^V) one can omit the exponentials in (j5l7p leading to 

V{r]) = K - / UW)di , K = V coth(^) . (14) 
Jo -^Co 

Note that K = corresponds toV* = (see ©) for small a. Now V can be replaced 

by the constant Vg in (U^V) (the rest is an odd function) leading to 

and from (jl3j) we infer 
2.1 Isolated pulses 

We consider now isolated pulses for which the wavelength L is large compared to the 
decay length of the inhibitor l/(ea). In Figures 1 and 2 we show a comparison of our 
results Eqs. ()l()pi2pi3j) for fo,w and cq, with a direct numerical simulation of Eqs. dlj. 
The pulse shape, the critical bifurcation point ec and the behavior near the saddle-node 
bifurcation of the amplitude /o(e) and of the velocity co(e) are very well recovered. Note 
the square-root behavior near the saddle node. 

Let us discuss some systematic features of the isolated pulses at criticality depending 
on the equation parameters a and Us, which can be extracted from our approach. Solutions 
/o of (jl(J|) exist when the discriminant — 4AC is positive. The amplitude at the saddle 
node fc is determined by the condition — AAC = 0. The corresponding bifurcation 
parameter ec can then subsequently be determined using Eqs. (fT^ and (fT^. Note that 
the saddle node, which occurs at the maximal bifurcation parameter ec, is not given by 
the relation B^ - AAC = since the condition of maximal e is different from that of 
maximal (see for example (fTT)|) ). In Fig. 3 we show ec and the corresponding amplitude 
at the saddle node fc as a function of Ug for a = (continuous line) for the isolated pulse 
with Vg = The points are results from a full solution of the 1-dimensional version 

of the ODEs ()2I3|) . In our approximation the limit of excitability (that is the maximum of 
Us for which 0) is = ^2(81 - 50/^2-9^81 - 50^2)750 = 0.4745, which is close 
to the exact stall-value = 0.5. Note that Uc is independent of a, as it should. Note also 
that pulse propagation in the neutrally stable case = is possible. Moreover pulses 
are supported even for negative Ug which we have checked numerically. We mention that 
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for this calculation it is crucial to take into account the difference between maximizing e 
and maximizing B. 

In Fig. 4 we show ec and fc as functions of a for = 0.1. We see that our approxi- 
mation (continuous line) reproduces all features of the ODEs (points). Fig. 4 reveals that 
around a = 1.3 the saddle- node bifurcation ceases to exist. The velocity cq approaches 
zero for these values. This correlates well with the fact that in the full system pulses be- 
come delocalized around a = 1.25, i.e. front and tail of a pulse separate creating a domain 
with u = 1, V = 1/a, which represents a locally stable stationary state of the system. 
As a matter of fact, the inverse pulse width w diverges here for the test function approach. 



2.2 Periodic wave trains 

Even if a given set of equation parameters allows for propagation of a single solitary pulse, 
the system may not necessarily support a wave train consisting of several of such pulses. 
As a matter of fact, if the distance L between two consecutive pulses of the train becomes 
too small, the pulses run into the refractory tail of the preceding pulse and consecutively 
decay. The critical wavelength Lc is a lower bound for the wavelength for the existence of 
periodic wave trains. On can also think of keeping L fixed and, as before, vary e. Then 
the saddle node €c{L) is a monotonically increasing function. 

One can calculate (or ec) essentially as before, except for the complication that due 
to Vs (fT^ . the expressions for /o,w and Cq, (|M j) .(fT ^ and ^7\i. cannot be cast in closed 
form depending only on G. This is true even in the limit of small a. However, given the 
equation parameters e,a and Ug one may obtain Lc numerically as a consistency relation 
requiring that at each L there exists a G = 1/ (cqw) so that the value for cq obtained by 
solving Eqs. (jlUll2p and using the relation cq = Q/w, matches the value for cq obtained 
by solving Eq. ()13|) . We obtain very good agreement between our test function approach 
and the numerically obtained values for the critical wavelength Lc. In Fig. 5 we show 
a comparison of the values obtained by integrating the full system and Q with the 
calculation of the test function approach as described above. The critical wavelength Lc 
diverges when e approaches ec where the saddle node of the localized pulse (i.e. L = oo) 
causes propagation failure of isolated pulses (see Section 2.1 and Fig. 2). 

In the remainder of this section we will discuss the limit of large values of L, i.e. small 
perturbations to the saddle node ec of the isolated pulse. This causes small shifts of the 
critical e, the amphtude /o and velocity Cq when compared to their respective values in 
the case of isolated pulses and L = oo. For simplicity, we restrict the calculation to the 
particular limit of small a. 

We write 

Vs = ^Vl+r , r = sinh-2(^) ^ 4exp {-aeL/co) (17) 
2 2co 

and expand in terms of small r. The correction to the isolated pulse r is connected to the 
exponential tail of the inhibitor from the previous pulse. It thus captures the interaction 
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between pulses mediated by the inhibitor. At leading order the shift with respect to the 
case L = oo of e^, /o and Cq at the saddle node must be proportional to r. For small a 
one finds at leading order 



(ec(oo) - ec(L))/ec(oo) 



r = 7exp (— aeL/co) , 7 = 4. 



(18) 



To see this note that the parameters /o, w, Cq are entirely given as a function of i? : = 
QVs via Eqs. (UHl) , (|I21) and Hence, if we multiply Eq. ^ by V^, we have four 

equations to determine the parameters fo,w,Co and R. We expand fo,R and ec with 
respect to r around the solutions of the isolated pulse. In particular, we write e^L) = 
ec{oo) + rei. The first-order correction ei can be entirely determined using 



where we deliberately ignored the subscripts to denote that / and c need to be expanded 
in r. The saddle-node condition de/df = implies that the derivative of the right-hand 
side of (I19j) with respect to / vanishes at r = 0, and we obtain the first-order correction 
61 leading to ()18|) . This result is confirmed by our numerical simulations. We have 
determined the shift of ec due to finite wave length L for a — > and = 0.1 numerically 
from the ODEs and find 7 = 4.3. This shows the accuracy of our test function approach 
for the saddle node shift. The leading-order approximation for the saddle node behavior 
is good down to about aL = 40. For a = 0.22 we find numerically 7 = 5.5. 

In Sec. |3]we will touch on some questions of stability of the wave train solutions near 
the saddle node. 

3 Growing velocity and retracting fingers in two di- 
mensional excitable media 

In this Section we develop a 2-dimensional extension of the test function approach. We 
study isolated finger solutions, i.e. solutions which in the —y direction go over into an 
isolated pulse (moving with velocity Cq in the x direction) and rapidly decay to zero in the 
+y direction. In the y coordinate they may be regarded as fronts, which grow or retract 
with velocity Cg (see Fig. 6 for a retracting case). We derive an explicit formula for the 
growing velocity Cg. We now investigate the full 2-dimensional system Q and Q in a 
frame moving with velocity (cq, Cg). 

We introduce a product ansatz for the activator field 



with test function U (77). This approximation neglects possible curvature at the tip. Again, 
we avoid any further uncontrolled approximation and solve for f{y) and v{x,y) in a 
systematic way. Note that f{y) replaces the constant /o in (jH). The solution of Eq. Q 
with the ansatz (|20|) can be written explicitly as 




(19) 



u{x,y) = f{y)U{ri) . 



(20) 



v{ri,y) 



e 



rn 
J 00 



e^^i^-^)f{h^{s - r/) + y/^))U{Us + {^?v))ds (21) 
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where A = 1 + {-^Y- Further calculations with this full expression appear prohibitive. 
Since we are mainly interested in the reversal point of Cg we resort to an (asymptotic) 
expansion in powers of c^, restricting ourselves here to the first two terms. From (|2ip one 
obtains 

y) = [9o{y)Vo{v) + cMy)yM + o{cl)) , (22) 

where 



9o{y) = f\y) and g,(y) = f\y). (23) 

Vo(?7) coincides with V{ri) of Eq. (0) with V^* = and L = oo (or, for small a, with Eq. 
Jill)) and 

Vi{r]) = ^e"®^ r e-'^®^Vo(^') d?]' . (24) 

cqw 

Note that the first-order correction of the inhibitor (j^^ can also be obtained by inserting 
the ansatz ()22p into Q and solving for successive orders of Cg. 

Repeating the procedure that led to Eq. (jH)) for /q, i.e. multiplying Eq. (j21) with U{r]) 
and integrating over r], we obtain using 



Df" + Cgf + F{y)f = (25) 

with 

F{y) = l-Dw^U'^) + {U'il - f{y)U)[f{y)U - Q{f{y)V, + c,/'(y)\/0 - u,])/{U^). (26) 

Note that for f{y) = /o±, where /o± are the solutions of the quadratic equation (fTUj) . we 
have F{y) = 0. 

We first neglect the higher-order correction of the inhibitor field Vi{y). Then F is a 
quadratic form in f{y) with zeros /o±, and Eq. ()25p can be solved exactly with the ansatz 

f'{y) = «/(/ - /o) , (27) 

which states that far away from the tip f{y) is constant and takes values or /q. The 
constant of proportionality a can be determined. Using ()27|) we obtain for the growing 
velocity 



2"V W) 



= \l-x r (/o+ - 2/o-) . (28) 



The point c^o = is fixed by the condition /o+ = 2/o_. The value for e where c^o = 
which we denote by e^, matches very well the value obtained by numerically integrating 
the full 2-dimensional system Q and Q. However, the behavior for nonzero growing 
velocities is not captured by (j28|) . In fact, the slope dcg^jde close to the reversal point is 
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too small by an order of magnitude for the parameters used in Fig. 1,2,5,6 when compared 
to the values obtained by the full 2-dimensional simulation. 

To obtain the correct slope we need to take into account the correction Vi{ri) of the 
inhibitor field. Including the first-order correction Vi(?7) in Eq. (j^^ we solve for Cg 
analogously to the determination of cq in Section 2 by multiplying Eq. by f'{y) and 
subsequently integrating over y. To 0{cg) we obtain 

_ 1 

" '''l + lGofo + ^,GJ^ ' ^^^^ 

where 

^0 = -©^^ and Gi = e^^. (30) 

As expected the higher-order corrections Go and Gi do not change the value of e^, but 
change the slope of dcg/de. In Fig. 7 we show a comparison of the test function approach 
and of Eq. (j^^ with numerically obtained data. The correspondence close to = is 
striking. To obtain better agreement further away from Cg = one would have to include 
higher-order terms in the expression ()22j] . 



4 Summary and discussion 

We have developed a non-perturbative method to study critical wave propagation of 
single pulses and periodic wavetrains in 1 and 2 dimensions. The method is based on 
the observation that near the bifurcation point the pulse shape is close to a symmetric 
bell-shaped function. A test function approximation, optimizing the two free parameters 
of a bell-shaped function, i.e. amplitude and width, allows to calculate the wave speed of 
a critical and close-to-critical pulse. We were able to study propagation failure of isolated 
pulses and wave trains, and moreover the test is capable of capturing more general features 
such as the transition from excitability to bistability. We have also performed our test 
function method with a more general class of nonsymmetric test functions to calculate the 
pulse parameters and velocities. It turns out that near the saddle node the asymmetry 
indeed becomes irrelevant. 

We extended our method to two dimensional situations, and used it to study broken 
fronts. Depending on parameters these fingers may either retract or sprout and start 
spiraling. We studied the growing velocity of a critical finger whose growing velocity 
is close to zero. Our test function method combined with a separation ansatz for the 
two-dimensional finger tip yields analytical expressions for the growing velocity which 
depend only on the equation parameters and which are for large parameter ranges in 
good agreement with the numerically obtained values. We elaborated on the importance 
of the inhibitor field for the growing velocity. 

Let us finally mention an interesting observation for one-dimensional wave trains. Be- 
low some (rather large) pulse separation (wavelength) the wave train looses stability not 
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via a saddle node but via a Hopf bifurcation. This may be seen in numerical simula- 
tions of the full system (0) and Q (and also of other excitable media equations such 
as the Fitzhugh-Nagumo equation) and in analytical approximations, either based on a 
time-dependent generalization of the test function approach or on a systematic reduction 
scheme valid near the saddle node. At the codimension 2 point where the saddle node 
of the wave train coincides with the Hopf-bifurcation the Hopf frequency becomes infi- 
nite. This time-dependent phenomenon can of course not be captured within the current 
stationary algebraic framework. This phenomenon and its implications will be published 
separately. For completeness we have included a time-dependent extension of the test 
function approach which allows for Hopf bifurcation in an Appendix, this seems to be 
related to the phenomenon of alternans [7] which has been studied in a different parameter 
regime far away from the saddle node jH]. 

Another interesting problem we plan to address is the selection of the wavelength (or 
pitch) of spirals. In simulations we found that near e^, where growth of fingers becomes 
small, the selected wavelength diverges. Kinematic theory HH 1121 addresses this 
problem. Kinematic theory provides, in principle, a relationship between the rotation 
frequency of a spiral and its core radius. However, it fails to provide expressions for either 
of the two which only depend on the equation parameters. A connection of our theory 
with kinematic theory is planned for further research to fill this gap. 
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A Appendix: Time-dependent system for wave prop- 
agation of pulses and periodic wave trains in one- 
dimensional excitable media 

In this appendix we present the time-dependent calculation for our test function approach. 
For simplicity we restrict ourselves to small values of a. This means that the width of 
the activator pulse u is small compared to the width of the decaying inhibitor field v. We 
include now temporal dependency of the pulse variables /o, Cq and w. We study pulse 
trains and note that the localized pulse can be obtained in the limit L ^ oo. We hence 
choose u of the form 

N 

u{x) = J2fnm{Vn) . (31) 


The sum extends over all N pulses. We defined rjn = Wn(t){x — (pnit)) and U{r]) is chosen 

2 

as a Gaussian e"'' as above. We allow for individual dynamics of the pulses characterized 
by the amplitude fn(t), the inverse width Wn(t), and the position (pnif)- For a stationary 
wavetrain the and Wn will be constant and all equal, and 0^ = n{c/L)t, where c is the 
velocity. We will restrict ourselves to the situation where the perturbations around such a 
state are small and slowly varying. We insert the ansatz (jHTj) into the general expression 

v{t,x) = -€ [ e-"'^''-'^u{t\x) ,dt' 



obtained from the second Eq. (QJ). Assuming that the (temporal) inverse pulse width 
{wnCn)~^ is small compared to the inhibitor decay time (ea)~^, we obtain the following 
expression valid in the vicinity of the n = pulse, which is assumed to pass the origin at 
t = 

v(r^,) = eJ^[l-{U) - rrfr/W)]+e(f/)Ve"'-^/^'^ + So /" rfr/V WI32) 
wqCq 2 Jq wici Jo 

where the brackets indicate integration over the whole //-domain, and 

So = -e{dt[fo/{coWo)])/{wofo)- (33) 

In order to determine /„, Wn, and c„ we again project Eq. (0) onto U{ri), riU'{r]), and 
U'{ri). After combining the first two equations appropriately one obtains 



- ^) = -ns{U') + |(1 + ns){U')fo - l{U')fo 
fo Wo 6 4 



-{{U') - l{U')h)V, - (^(t/3) - ^(f/^)/o)| , (34) 
{U')^ = 2Dwl{U'') + \{l+u.){U')fo - l{U')f^ 

Wo O 

+lmfoK - iliU') - ^(f/')/o)? , (35) 
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where 



1 /o 



oo 



fl 



2coWo 



WlCi 



(36) 



1=1 



The third projection gives an algebraic relation 



CqWq Z 



(37) 



These equations are written for the pulse n = 0. They apply correspondingly to the other 
pulses. 

First consider the case of a stationary, isolated pulse {L = oo) where Sq = and the 
sum in Eq. ()36|) vanishes. The resulting algebraic relations are easily solved, and corre- 
spond in the limit of a ^ to Eqs. (fTUI) . (fT^ . (fT^ . With reasonable initial conditions, 
simulation of the ODEs (jH^ . (jH^j) . together with (jHTf) leads to the same result as described 
in Section 2. The ODEs (j34l35j) relax to the stationary values obtained in Section 2. 

This time- dependent test function approach allows to go beyond the stationary bi- 
furcations discussed in Section HI As mentioned in Section 2 the saddle-node bifurcation 
related to propagation failure for well separated pulses transforms into a subcritical Hopf 
bifurcation via a Takens-Bogdanov point when the pulse separation is reduced below a 
critical value Lc- This work will be published elsewhere. The Hopf-bifurcation can be 
captured within the ODE-system (jSH), and and some preliminary simulations 
have been done. 
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Figure 1: Activator field u of a steady front witli e ~ ec at e = 0.0485. The points depict 
the solution obtained by numerically integrating Eqs. (0). The continuous line is the 
theoretical curve obtained with the test function approach. Note that here e is not even 
very close to ec ~ 0.049. 
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Fig. 2 
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Figure 2: (a): Amplitude fo as a function of e. (b): Velocity cq as a function of e. The 
parameters used were D = 3.0, a = 0.22 and Ug = 0.1. The crosses depict the numerically 
obtained values of integrating the full system Q and Q, the lines depict the stable and 
unstable branch calculated with the test function approach, i.e. by using Eqs. ()10I12I13|) . 
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0.4 




Figure 3: Behavior at the saddle node of an isolated pulse for a = 0. (a): The critical fc 
versus Ug- (b): Critical amplitude ec versus Us- Points correspond to simulations of the 
full 1-dimensional version of (j2l3|) . The continuous line corresponds to the test function 
approach. 
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Figure 4: Behavior at the saddle node of an isolated pulse for Ug = 0.1. (a): The critical 
ec versus a. (b): Critical amplitude fc versus a. Points correspond to simulations of the 
full 1-dimensional version of (j2l3|) . The continuous line corresponds to the test function 
approach. 
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Fig. 5 
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Figure 5: (a) The critical wavelength Lc as a function of e. Crosses depict the values 
obtained by direct simulation of the full system Q and ©• Squares depict the values 
obtained by the test function approach described in Section 3. Parameters are again 
D = 3.0, a = 0.22 and = 0.1. In (b) the same numerical results as in (a) are presented 
but here the data-points corresponding to the test-function approach are shifted along 
the e-axis such that the saddle nodes at Lc = oo (see Fig. 2), coincide. 
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Fig. 6 
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Figure 6: Contour plot of the activator m of a retracting finger at different times. 
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Figure 7: Growing velocity Cg as a function of e. The crosses depict the values obtained by 
direct numerical integration of Eqs. (j21) and Q- The continuous line shows the theoretical 
curve using the test function approach of Section 4.1. We have shifted the curve along 
the e-axis by the difference of the e-values A e = 0.001972 for the saddle nodes obtained 
by the numerical simulations of the full system (j2)) and Q, and the test function approach 
of Section 2. 
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